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ABSTRACT 

In this work, several algorithms based on higher-order moment (autocorrelation) 
matching of single hydrophone element data have been developed and tested on real 
transient data sets. Of particular interest is the success and robustness of the Frequency- 
domain Autocorrelation Matching (FACM) algorithms in the presence of environmental 
mismatch, signal mismatch, and noise, for different signals in an unknown environment. 
Recently acquired data was analyzed for signal variability in terms of spatial coherence 
of phones, beams, and modal structure. The ability to localize using these higher-order 
moment matching algorithms was compared to the spatial structure of the signal, the 
placement of the receiving elements, and the signal variability. 

This work suggests that the FACM algorithms are strongly dependent on the 
source-receiver relative positions, and on the uniqueness of the signal vertical structure. 
It is also shown that their performance increases with the number of multipath arrivals 
and, therefore, with the range. More importantly, the localization results obtained with 
raw linear frequency modulated (LFM) signals seemed to be as useful as the ones 



DUDLEY KNOX LIBRARY 

naval postgraduate school 

tevNfEREY Ca 93943-5101 



obtained from matched-filtered data. 



TABLE OF CONTENTS 



l. INTRODUCTION 1 

A. THE NUWC-TSP EXPERIMENT 2 

B . THESIS OBJECTIVE APPROACH AND OUTLINE 7 

II. THE MONTEREY -MIAMI PARABOLIC EQUATION PROPAGATION 

MODEL 9 

m. THE FREQUENCY-DOMAIN AUTOCORRECTION MATCHING 

(FACM) LOCALIZATION ALGORITHM 17 

IV. SIGNAL PROCESSING AND ENVIRONMENTAL ASSESSMENTS 25 

A. MATCHED-FILTER PROCESSING 25 

B. PLANE- WAVE BEAMFORMING 27 

C. MODAL ANALYSIS 33 

D. SPATIAL COHERENCE (CORRELATION) 37 

V. LOCALIZATION RESULTS 41 

VI. CONCLUSIONS 49 

A. SUMMARY OF FINDINGS 49 

B. RECOMMENDATIONS FOR FUTURE WORK 50 

LIST OF REFERENCES 53 

INITIAL DISTRIBUTION LIST 55 



vii 













viii 




LIST OF FIGURES 



1 . 1 The NUW C-TSP Experiment Geographic Location 3 

1 .2. The VLA Configuration 4 

1 .3. The NUWC-TSP Experiment Configuration 4 

1.4. The Transmitted LFM Burst (Left) and the Received Signal at VLA 

(Right) in the Time-Domain. Note That the Axis Are Different on Each Plot.... 5 

1.5. Frequency Spectra of Transmitted Chirp (Left) and Received Burst (Right) 6 

1.6. Sound Speed Profiles (SSP’s) Collected During Experiment 7 

4.1. Examples of a Raw (Left) and an MF Burst (Right) in the Time -Domain 26 

4.2. The Frequency Spectra of a Raw Burst (Left) and of an MF Burst (Right) 26 

4.3. Beamformer Outputs for Scenario 1 1 (Range 460 m, Depth 45.7 m) 30 

4.4. Beamformer Outputs for Scenario 16 (Range 2194 m. Depth 27.4 m) 31 

4.5. Beamformer Outputs for Scenario 9 (Range 4650 m, Depth 45.7 m) 32 

4.6. The First 10 Normal Modes (f = 900 Hz). The Bottom Depth is 89.9 m 34 

4.7. Modal Amplitudes vs. Reduced Time-Seen. 1 1 (Range 460m, Depth 45.7m)... 35 

4.8. Modal Amplitudes vs. Reduced Time-Seen. 16 (Range 2194m, Depth 27.4m). 35 

4.9. Modal Amplitudes vs. Reduced Time-Scen.9 (Range 4650m, Depth 45.7m)... 36 

4.10. Phone-to-Phone Peak Correlation - Scenario 1 1 (Range 460 m) 38 

4.1 1. Phone-to-Phone Peak Correlation - Scenario 16 (Range 2194 m) 39 

4. 12. Phone-to-Phone Peak Correlation - Scenario 9 (Range 4650 m) 40 

5.1 Ambiguity Surfaces for Scenario 9 - (a) with Raw Data from Hyd. 16 
(Depth 54.8), and (b) with Raw Data from Hyd. 2 (Depth 76.3) 

The White Dot Depicts the Actual Source Location 44 

5.2 Ambiguity Surfaces for Scenario 9 - (a) with Raw Data and (b) with MF 
Data from Hyd. 16 (Depth 54.8 m). The White Dot Depicts the Actual 

Source Location 45 

5.3 Ambiguity Surfaces for Scenario 16 - (a) with Raw Data and (b) with MF 
Data from Hyd. 16 (Depth 54.8 m). The White Dot Depicts the Actual 

Source Location 46 

5.4 Ambiguity Surfaces for Scenario 1 1 - (a) with Raw Data and (b) with MF 
Data from Hyd. 16 (Depth 54.8 m). The White Dot Depicts the Actual 

Source Location 47 



IX 





X 



I. INTRODUCTION 



In the past, passive sonar systems were only able to provide bearing information, and 
the estimation of source range relied on analytical, over-simplified techniques. Starting in 
the mid-60’s, however, with relatively advanced computational resources, the received 
signals have been processed in order to also provide information about the distance to the 
source. Today, one of the most used and successful research techniques is matched-field 
processing (MFP) where a comparison is made between the parameters of the received 
signal and those of signals generated by a synthetic source virtually positioned at each point 
in the search grid. A match is obtained when the correlation of the synthetic and the true 
signals are at a maximum. 

Transient signals are especially important for naval purposes, since they can represent 
several instances of a ship’s routine, such as firing a torpedo, opening of hatches, starting 
of pumps, etc. Therefore, the Transient Localization Project at the Naval Postgraduate 
School has been studying, since 1993, several algorithms for localizing transient sources 
based on methods similar to MFP. 

In 1995, three distinct schemes were studied: (a) the fully coherent localization, which 
was considered impractical as it required an extremely precise (and unavailable) 
propagation model, along with accurate environmental descriptions; (b) the semi-coherent 
localization, which correlated only the slower varying amplitude terms of the signal, 
trying to surpass the limitations of method (a); and (c) the autocorrelation matching, which 
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compared the autocorrelations of the received and modeled signals using a simple Bartlett 
technique. The NPS results showed that method (b) was not successful in localizing the 
signals of naval interest, and method (c), while not perfectly consistent, was the most 
useful technique. According to Miller et al. (1996), autocorrelation matching has an 
advantage over the other methods, since its processing allows the notching of the noise 
energy and the possibility of a multi-transient gain. 

A. THE NUWC-TSP EXPERIMENT 

The transient signal processing (TSP) experiment was a field study sponsored by 
NUWC, and conducted both by NUWC and C&M Tech., Inc. on 28-31 July 1997. Its 
purpose was to provide data to study the influence of a highly variable shallow water 
environment on acoustic propagation, and the effectiveness of localization algorithms. 

The study area is shown in Figure 1.1 and as it is influenced by oceanographic and 
geological processes related to the shelfbreak, it provides a complex coastal environment 
for sound propagation. The continental shelf water is cold, and it presents large variations 
with location and season. The slope water, on the other hand, is warmer and influenced by 
the Gulf Stream.(Miller, C. W„ 1998) 
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Figure 1.1. The NUWC-TSP Experiment Geographic Location. 

Synthetic signals were generated and transmitted in a variety of scenarios which were 
differentiated by the source depth, source-receiver range, and bearing. Two projectors 
(ITC-104 and HX-188), a 32-element vertical- line hydrophone array (VLA), and data- 
recording equipment constituted the main assets of the experiment. The VLA 
configuration is depicted in Figure 1.2, and the experiment configuration for one of the 
scenarios is shown in Figure 1.3. 
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Figure 1.2. The VLA Configuration. 
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Figure 1 .3. The NUWC-TSP Experiment Configuration. 
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It is important to note that hydrophones 17 and 18 were non-operative throughout the 
experiment. In this work, we were only interested in the linear frequency modulated 
(LFM), short duration (100 ms) signals - the transmitted chirp and a sample of a received 
burst can be seen in Figure 1.4. Their frequency spectra are observed in Figure 1.5. 

The processed data corresponded to three particular cases: scenario 9 (source located 
at 4650 m range, 45.7 m depth), scenario 11 (source at 460 m range, 45.7 m depth), and 
scenario 16 (source at 2194 m range, 27.4 m depth). The information on source location 
was provided by NUWC, based on ship’s log and GPS data. While some uncertainties are 
present in this information, inherited from the processes used to gather it, we assumed it as 
our reference. The recorded tapes also contain array tilt information, which was neglected 
in this work for simplicity. 




Figure 1.4. The Transmitted LFM Burst (Left) and the Received Signal at VLA (Right) in 
the Time-Domain. Note That the Axis Are Different on Each Plot 



5 





Figure 1.5. Frequency Spectra of Transmitted Chirp (Left) and Received Burst (Right). 

In order to support the analysis of the experimental data, six conductivity-temperature- 
depth (CTD) casts were made. Two of them did not present useful data and were 
disregarded. The other four provided the sound speed profiles (SSP) shown in Figure 1.6, 
which also displays the date and time (Zulu) they were taken. For processing, however, we 
used only SVj97_08 (applied to scenario 9), and SVj97_12 (applied to scenarios 11 and 
16), because they were gathered at times closer to those of the studied runs. Still, it is 
important to note the large variations, particularly near the lower part of the thermocline 
' between 20-35 m. The associated pycnocline supports the propagation of nonlinear soliton 
waves which travel from the edge of the slope onto the shelf. Such oceanographic features 
can create localized regions of high sound speed contrast as strong as 25 m/s over a 100 m 
horizontal range in the direction of the soliton propagation (Chiu et al., 1997). 
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Figure 1.6. Sound Speed Profiles (SSP’s) Collected During Experiment. 

B. THESIS OBJECTIVE APPROACH AND OUTLINE 

The objective of this thesis is to examine the influence of shallow water acoustic 
variability on frequency-domain autocorrelation matching (FACM) localization. In order 
to perform such analysis, we processed the acoustic VLA data, identified the LFM bursts, 
and high-pass filtered them to reduce the array strum influence (cut-off at 50 Hz). The 
resulting bursts, as well as the replicas generated by the MMPE propagation model, were 
applied to a Bartlett-type processor based on matching broadband signal autocorrelation 
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localization algorithm). The acoustics variability study followed with the plane-wave 
beamforming, modal amplitude structure, and spatial coherence analyses. 

The remainder of this thesis consists of five chapters. Chapter II describes the 
Monterey-Miami Parabolic Equation (MMPE) propagation model. Chapter IE discusses 
the Frequency-domain Autocorrelation Matching (FACM) localization algorithm. Chapter 
IV describes the techniques used to evaluate the environmental influences and process the 
signals. Chapter V presents the localization results. Chapter VI concludes this study with a 
summary of findings and recommendations for future work. 
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II. THE MONTEREY-MIAMI PARABOLIC EQUATION 
PROPAGATION MODEL 



The Monterey-Miami Parabolic equation (MMPE) acoustic propagation model 
(Smith, K., 1996) was used to generate the replicas which simulated the received signals 
from a synthetic source on the gridpoints of a 2-D search space (depth x range). It allowed 
us to predict the arrival structure at the VLA due to a transient- like, broadband, point 
source. In this chapter, we introduce the general theory behind the parabolic equation 
model, as well as the method used for its implementation - the Split-Step Fourier (SSF) 
method (Hardin and Tappert, 1973). 

The parabolic equation method for computing underwater acoustic propagation was 
introduced by Tappert (1977). To begin the parabolic approximation development, 
consider the representation of the time-harmonic acoustic field in cylindrical coordinates 



As the model is based on an approximation to the Helmholtz wave equation, substituting 
Eq. 1 into the latter in cylindrical coordinates leads to (Jensen et al., 1994) 



P(r,z,ty,f,t) = pj{r,z,§)e l2nfi 



( 1 ) 



2 



2 




( 2 ) 
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where n(r, z, (t>)= — — is the acoustic index of refraction, c„ is the reference sound 

c(r, z, <j>) 



CO 



speed, c(r, z, <{>) is the acoustic sound speed, and k = — is the reference wave number. 

u C 
u o 

In this derivation, however, we neglected density variations, which could be 
incorporated into a new index of refraction without any loss of generality. Note also that 
the environmental characteristics are within c(r, z, <|>), and the reference pressure level 
P 0 is defined as the pressure amplitude at a distance R 0 = 1 . 

By assuming the ocean acts as a waveguide with a cylindrical coordinate system, 
acoustic energy is primarily propagated outward from a source in the horizontal direction. 
The pressure field, therefore, can be approximated by 



where the slowly varying envelope V|/y(r, z, <j>) modulates the outgoing zero-th order 
Hankel function of first kind H^ iy >{k 0 r ) . 

In the far field (k 0 r » 1 ) , we can approximate the Hankel function by the first term of 
its asymptotic expansion (Gradshteyn and Ryzhik, 1994), 



P/(r, z, <(>) = \| fj(r,z,^)H^(k 0 r), 



( 3 ) 




( 4 ) 
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and Eq. 3 becomes 



Pfr, z, <]>) = P 0 Jyyjir,z,<l>). (5) 

Note that Eq. 5 is normalized such that at r = R Q , |vyj = 1 and (pyj = P Q . It expresses 
the relationship between the PE field function z, <|>) and the acoustic pressure 
Pj(r, z, <t>) . By substituting it into Eq. 2 and dropping the source term on the right hand 
side, we obtain 



+ y *[&*-»*■&}' ■ °' <6) 



If we neglect the azimuthal coupling and the far field terms, Eq. 6 can be written as 



F57 - <^ + £ Vv. 



(7) 



where 



T =-M— 
op 2k 0 {.dz 2 



( 8 ) 



and 






(9) 



2 

Note that this is only valid for « \j/ } once we used the uncoupled azimuth 



approximation. 
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Eq. 7 is known as the standard parabolic equation (SPE), and the accurate solutions are 
limited to a half beam width of 10° to 20° for the propagation angle. However, in order to 
extend this limit to 40° ~ 70° , a higher order wide-angle parabolic equation (WAPE) 
approximation (Thomson and Chapman, 1983) can be used. It also features less sensitivity 
to the choice of k 0 , and to the phase errors in typical deep ocean conditions, relative to the 

SPE approximation.(Jensen et al., 1994 and Chin-Bing et al., 1993) 

The WAPE operators are defined as 

r = 1 d 2 

TwAPE ~~k^ 

and 

U\vape = ~( n - 1 ) • ( 11 ) 

In order to numerically solve the parabolic equation, the MMPE uses the split-step Fourier 
(SSF) method. This algorithm integrates the solution in range by applying the operator 
U wape i* 1 t ^ le z-domain., and the operator T WAPE in the -domain. The latter is defined, 

in the k z -domain, as 



1 + 



k 2 dz 2 



+ 1 



( 10 ) 



TwAPE(k 7 ) = 1 - 



k}') 



1 -"3 

k 2 

K oJ 



( 12 ) 



Note that both operators are just scalar multipliers, and may be applied independently. 
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The approximate solution for \\f is iterated in range according to the following 



expression 



\|/(r + A r, z) = e 



-ik 0 Ar U wAPE ( r » £) 




- ik 0 ArTwAPE(r , k z ) 



[FF7Xy*(r,z))]*), (13) 



where the FFT calculation is consistent with the form described in (Press et al., 1988). A 
forward FFT and double conjugation emulate the inverse FFT. The model outputs the field 
functions \\i (both magnitude and phase), referenced to a unit magnitude at r = lm, and 
computed at the spatial grid points. 

The source assumed in our model had a 600 Hz bandwidth centered at 900 Hz. 
Therefore, the acoustic field was computed for each one of the discrete frequencies 
(N=512) of the considered BW, and the model properly represented the broadband 
acoustic field. 

As the single component of the general field can be expressed as Eq. 5, and assuming 
a windowed (Hanning), normalized source amplitude, the time-domain complex pressure 
field can be written as 





(14) 



13 



ik r ilKf 7 

Note that the overall phase factor e ° = e ° can be neglected if we consider that the 



arrival times are given in terms of “reduced time” T 




, then 




(15) 



This means that the time domain is heterodyned around t = — , and the use of a reduced 

time does not have any influence on the autocorrelation function. The model also 
heterodynes the signal by shifting the center frequency to zero (d'.c.), in order to reduce the 
computational load associated with large transform sizes. This procedure does not 
introduce any consequence to our algorithm. 



The MMPE model, a FORTRAN® code, is the latest version of what was formerly 
known as the University of Miami Parabolic Equation (UMPE) model (Smith and Tappert, 
1993). It requires a series of input files with strict formats, which properly represent: 

• environmental data (sound-speed profile, bathymetry of water/bottom interface, 
bathymetry of the deep layer beneath the water/sediment interface, acoustic 
parameters of this “deep” layer, and acoustic parameters of the medium just below 
the water/bottom interface); 

• source data, with two types available, wide-angle and vertical line array (steering 
allowed). We only used the wide-angle, which approximates the point source. The 
center frequency, bandwidth, and number of discrete frequencies (a power of two, 
because of the FFT) are also required. 

All these input files are specified in a main one, the “pefile's.inp”, where several other data 
(name of output files, grid size and grid steps, c 0 , size of vertical FFT, etc.) are contained. 
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The output is presented in a single binary file composed of a header (with the information 
needed for post-processing), and the complex PE field function y at selected grid points 
for every discrete frequency. 
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III. THE FREQUENCY-DOMAIN AUTOCORRECTION MATCHING 
(FACM) LOCALIZATION ALGORITHM 



The algorithms studied in the Transient Localization Project at Naval Postgraduate 
School are designed for data from a single receiver hydrophone and a point source. The 
basic concepts of the correlation functions are outlined in Bendat and Piersol (1971), and 
the routines are described in Miller et al. (1996), Pierce (1996), and Hager (1997). From 
Brune (1998) we have that: 



Localization algorithms may be considered generalized 
beamformers in which the plane wave replicas have been replaced by 
more complicated replicas of the acoustic propagation (e.g., modes, 
beams, or the vertical pressure field). The algorithms, usually referred to 
as processors, are in most cases based on a Hermitian quadratic product. 
The exact form is determined by the constraints that are put on the 
processor output. 



As discussed in Chapter II, replicas are simulated received signals of synthetic sources 
located at the points of a search grid. Generated by a propagation model, they are matched 
with the received signal. Considering that the source functions are equal and a perfect 
model is used, the exact match should occur where real and synthetic source locations 
coincide. An ambiguity surface usually represents the localization algorithm results. 

The reciprocity principle states that, in a time-invariant and homogeneous 
environment, the acoustic pressure at location A due to an omnidirectional source located 
in B is equivalent to that of location B due to an omnidirectional source located in A. 
Using this principle, we reduce the computational load required to obtain localization. 
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since we are able to set the propagation model to generate replicas from the receiver 
position to all source possible locations in the grid. 

In this chapter, we discuss the derivation of the frequency-domain autocorrelation 
matching (FACM) localization algorithm. Autocorrelation matching was originally 
designed for use in the time-domain (TACM), therefore it is more convenient to start its 
study treating the problem in time, and then to convert the results to the frequency domain. 
Previous investigations done by de Kooter (1997) showed that the TACM algorithms 
present very small footprints and large phase errors at higher frequencies, being useful for 
very low frequency signals. It was also observed that matching the time-domain 
amplitude-squared signals, which seemed relatively similar, could be more effective and 
robust. Therefore, as we will describe in the following paragraphs, an almost natural 
further development led to the frequency-domain algorithm. We begin with a complex 
pressure time series P(t), which represents a detected transient arrival. 

The signal autocorrelation in the time-domain can be expressed by 




(16) 



or, in terms of the frequency-domain response, by 




(17) 



where 




(18) 
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Then, the normalized autocorrelation becomes 



Cppi't) — 



Tpp('t) 

f M 2 df 



A similar derivation applied to the predicted replica signal R(t) leads to 



( 19 ) 



c ss (t) = ■ ( 20 ) 

} k/il df 

— oo 

The autocorrelation matching algorithm is based upon the inner product of these two 
quantities and, normalized, in time-domain, it can be expressed by 



J Cftft(x)Cpp(x)dx 

— oo 

[f |C ffs (t)| 2 <iTj“ \C PF (T;fd% 



By analogy, in the frequency-domain it becomes 



J C ~RRW C ~PP^ d< b 

— oo 

[f \wm 



(21) 



( 22 ) 



where Cy p and Cy R are the autocorrelations in frequency-domain of the transient signal 
and of the replica, respectively. 
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Since 



f cl- R ($)Cy p W)di J)= f ^l^(0| 2 e'' 2 ^ r irj” dt'\P(t')\ 2 e i2K¥ ] 

— oo — oo I— — oo — i 1— — oo J 

= f dt\P(n\ 2 f dt\R( t f\f d<$>e i2KHt '- t) '} (23) 

— CO — OO 1— — CO —I 



= f \R(tf\P(tfdt 



Eq. 22 can be rewritten as 



J |*( 0 lV( 0 l 2 * 



Arp = 



[) \R(t)\ 4 dtf \P(t)\ 4 dt 

L oo CO — 



-,1/2 ■ 



(24) 



In order to reduce the influence of noise, we shall remove the zero-lag component of 
the autocorrelation function, which is equivalent to removing the mean of the squared 
amplitude of the time-domain signal. Therefore, we define 



|*'(t)| 2 = |tf(t)| 2 -<|/?(t)| 2 >, 



and 



(25) 



| P '( 0 I 2 ’= | P ( 0 l 2 -< LP ( 0 lV 

where < > t is an average over all times. Eq. 24 then becomes 



(26) 



f \R\tf\P\t)\ 2 dt 



Arp = 



1/2 



r— CO CO — l j 

|J \R\t)\ 4 dt\ \P\t)\ 4 dt 

1— — OO — CO — 

which is focused on the influence of matching the non-zero lag components. 



(27) 
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However, it still presents an ambiguity in (absolute) time which was not present in 



TACM, and we shall introduce a search parameter x to allow us to slide the replica signal 
in time, searching for the optimal match. The FACM function then becomes 



Arp - maxj 



f \R'(t + x)\ 2 \P'(tfdt 



[f \R'(t)\ 4 dtf \P\t)\ 4 dt~\ 

L — oo — oo —1 



1/2 



(28) 



where max x is the maximum value of the function for any value of the search parameter 
x. 



A MATLAB ® implementation of the algorithm based on Eq. 28 was used to generate 
all the ambiguity surfaces shown in this work. The sliding- x operation was performed 

2 2 

according to the following development. Let g{t) = |P'(*)I and h(t) = |i?'(t)l > then 



g(t) = f e- i2nft G(f)df , 

— oo 



(29) 



and 



h{t) = f e~ ,2Kft H(f)df 

— oo 



A« + t) = f f . (30) 

OO — oo 
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Therefore, the numerator of Eq. 28 becomes 



A(t) = J h(t + x)g(t)dt. (31) 



If h(t) is a real function, 



CO 

A(x) = | h*{t + z)g{t)dt 



= 1 



“ CO — 

J e i2n/t H*(f)e i2nfz df 


" CO ~ 

J e~ i2nrt G(f')df' 


CO 





dt 



JJ 

//' L 



/. ilTlfl f -i2n(f'-f)t 



H*(f)G(f')e i2nfz J e 



dt 



df'df 



CO 

= J H* (f)G(f) e' 2nfz df , 



knowing that 



(32) 



J e ^’-fi'dt = 5(f _/). • (33) 

—CO 

Hence, the search over x was in fact computed by a multiplication in the frequency 
domain and an inverse FFT operation represented by the numerator A(x). Prior 
investigations (Brune, 1998 and de Kooter, 1997) had shown that the FACM algorithm is 
better suited for larger bandwidths and higher frequencies, when compared to its analog in 
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the time domain, the TACM algorithm. Also, FACM is more robust against environmental 



mismatch and presents better peak-to-sidelobe levels and larger footprints. 
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IV. SIGNAL PROCESSING AND ENVIRONMENTAL 

ASSESSMENTS 



A. MATCHED-FILTER PROCESSING 

Matched filtering (MF) is a process where the received signal r(t) is correlated against 
a replica of the transmitted signal s(t), in order to obtain a better signal-to-noise ratio 
(SNR) and reduce the effective temporal structure of the source down to an ideal, coherent 
pulse. It is also known. as a conjugate filter, since its frequency response is basically the 
conjugate of the (transmitted) signal spectrum (Tolstoy, A., 1993). In this work, matched- 
filtering was used first to identify the short duration bursts (0.1s) in the 2-min. long data 
files. Second, but more importantly, to process the bursts and obtain data emulating the 
existence of a coherent source, instead of the actual non-coherent one. This also provides a 
baseline for the localization results since the model-generated replicas are based on a 
coherent source. Therefore, this is equivalent to localizing a transient with known source 
function. 

The process was implemented digitally in M ATLAB ® , according to the expression 



where p{t) is the matched-filtered output in time-domain, S*(f) is the complex conjugate 
of the transmitted signal in the frequency domain, and R (f) is the received signal in the 
frequency domain. 




(34) 
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An example of the MF output in time domain can be visualized in Figure 4.1, and its 
frequency spectrum in Figure 4.2. The raw data is also shown and, by comparison, we can 
observe the significant noise reduction. 
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Figure 4.1. Examples of a Raw (Left) and an MF Burst (Right) in the Time-Domain. 
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Figure 4.2. The Frequency Spectra of a Raw Burst (Left) and of an MF Burst (Right). 
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B. PLANE- WAVE BEAMFORMING 



Plane-wave beamforming is one of the many techniques used to distinguish particular 
features of a received signal while reducing the influence of noise - which is usually 
generated near the surface, and tends to occur in high modes or high angles of propagation 
(Jensen et al., 1994). It consists of decomposing the signal into plane waves by summing 
the outputs of spatially distributed sensors - in this work, eight hydrophones for each of 
the three sub-arrays in the VLA. 

In underwater acoustic propagation, much of the information about the environment 
and the relative source/receiver location is contained in the vertical structure of the 
acoustic field. Therefore, the ability to resolve the vertical structure of the arrival at the 
receiver can lead to estimations of the corresponding source location. 

A beamformer is especially useful because it allows us to give preference to one 
direction of propagation over another, implementing a spatial filter. In this work, it was 
used to observe the dominant beams when the signal was arriving at the VLA, providing 
information for a comparison of the localization algorithm results by the distinct arrival 
angle structures. 

An excellent discussion about the beamforming process is presented by Defatta et al., 
(1988) and its summary follows. Consider an infinite, homogeneous medium containing a 
line array of equally spaced elements positioned along the y-axis, and an undetermined 
number of remote sources. The output of an element located at the origin of coordinates 
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due to the Z 111 source is defined as S[(t ) . Assuming plane- wave propagation in the far- field. 



a source from direction 0 ; outputs 



( nJsin(0,)\ 

«.(») = s(»+ — 



( 35 ) 



where n is the array element index, d is the spacing of the elements, 0^ is the angle of 



arrival from the I th source, and c is the speed of propagation. To accomplish beamforming. 



we must apply weights 0 and time delays to the individual element signals, and then 
coherently sum them. Therefore, the output of the beamformer in the direction Q m is 



N- 1 

b m {t) = Z WjSA t + 

n = 0 



nd[s in(0 ; )- sin(0 m )]' 



( 36 ) 



To derive the directional response characteristics of the beamformer, we 'shall consider a 
continuous, complex plane-wave signal propagating across the array, 



I(0,f 



c ,u v 
S[ = e 



( 37 ) 



Due to the arrival angle 0^ , the sensor outputs 



i(£),t ink,dsm(Qi) 

e n (t) = e e 



( 38 ) 



2k 

where k, = = — is the wave number of the I th source. 

1 h c 
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The beamformer output is then 

N- 1 



, ,,, *'«¥ v ‘"Ml sin (9()-sin(0 m )] ( 0 ), « 

b m {t)= e Zj w n e = W(Q t )e , 

n = 0 



where 



N- 1 

W,) = Z 

n = 0 



ink l d[sin(Q l ) - sin(0 m )] 



(39) 



(40) 



is the spatial Fourier transform of the array weights. 

The beamforming operation applied in this work was implemented using FFT in 



MATLAB®. Single frequencies considered, the time-delay operations described above 
are equivalent to phase shifts that are linear functions of the element index. FFT 
beamforming is the implementation of these phase shifts through FFT operations. The 
processed results for all three sub-arrays, and the three distinct scenarios can be seen in 
Figures 4.3, 4.4, and 4.5. They describe the physical angle of arrival as a function of the 
reduced time, and provide an easy visualization of the multipath nature of the burst 
propagation. The results for sub-array 2, however, are noticeably affected by the lack of 
data from hydrophones 17 and 18. We can verify in Figure 4.3 (scenario 11, range -500 
m) that the energy is not as concentrated in the lower angles as in the other two figures 
corresponding to larger ranges. However, in all of them we observe the early arrivals of 
the lower angles, the short period in which the energy is contained, and the distinct 
multipath structure. 
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Tran am oston Loea (dB >o 1m) - Scenario 11, Sub-array 3 
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(a) Sub-array 3 



Tran am aa*on Loae (dB ra 1m) - Scenario 1 1 , Sub-array 2 




Tima (eoo) 



(b) Sub-array 2 




0 46 0 46 



(c) Sub-array 1 

Figure 4.3. Beamformer Outputs for Scenario 1 1 (Range 460 m, Depth 45.7 m). 
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(a) Sub-array 3 



TranerrtaaJon Loe» (dB re 1m) - Scenario 18, 9ub-array 2 




(b) Sub-array 2 
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(c) Sub-array 1 



Figure 4.4. Beamformer Outputs for Scenario 16 (Range 2194 m. Depth 27.4 m). 
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Transmission Lom (dB re 1m) - Scenario 9, 9uO— array 3 




Time (sec) 



(a) Sub-array 3 



Transmission Loee (dB re 1m) - Scenario 9. Sub-army 2 




Time (eao) 



(b) Sub-array 2 



Transmission Loss (dB re 1m) - Scenario 9, Sub-array 1 




Tlmo (aoc) 



(c) Sub-array 1 

Figure 4.5. Beamformer Outputs for Scenario 9 (Range 4650 m. Depth 45.7 m). 
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C. MODAL ANALYSIS 



Normal modes are another physical, frequency- dependent representation of the 
underwater acoustic wavefield. For a given sound channel, the acoustic modes represent a 
unique, natural set of standing waves oscillating along the depth axis (Chiu and Ehret, 
1994). 

The method involves solving a depth-dependent equation, and the complete acoustic 
field is then constructed by summing up contributions of each of the modes weighted in 
accordance to the source depth.The modal equation 



p<z) £l 



C 1 <W m (z) -I 


r 2 i 

0) 2 


Lp(z) dz J + 


2 K rm 

-C (z) 



0 . 



(41) 



with boundary conditions 4 / (0)= 0 (pressure-release surface at z=0) , and 



dz 



= 0 



z = D 



(perfectly rigid bottom at z = D) is a classical Sturm-Liouville eigenvalue problem, whose 



properties are well known (if we assume that p(z) and c(z ) are real). Some are: 

• the modal equation has an infinite number of solutions which are like the modes of 
a vibrating string; 

• the modes are characterized by a mode shape function l F m (z) (an 
eigenfunction ) and a horizontal propagation constant krm (eigenvalues). These 
constants are all distinct and analogous to a frequency of vibration; and 

• the modes are unique (orthogonal) and form a complete set, which means that we 
can represent an arbitrary function as a sum of the normal modes. 
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Based upon these and through a modal decomposition implemented in MATLAB, we 
were able to quantify the amplitudes of the first 30 modes as functions of the signal 
(discrete) frequencies and, by IFFT, of the reduced time. The SSP’s used were the ones 
discussed in Chapter I and depicted in Figure 1.6, applied to their respective cases. A 
visualization of the first ten modes of the assumed waveguide for a frequency f = 900 Hz 
is presented in Figure 4.6. The final results of the modal decomposition of the received 
signal for each studied scenario are shown in Figures 4.7, 4.8 and 4.9. Note the energy 
concentration in the lower modes, especially scenarios 9 (5000 m) and 16 (2500 m). This 
is consistent with stripping of higher modes at longer ranges. The energy observed in the 
higher modes must be treated with skepticism since the array does not sample those modes 
very well. The lack of tilt information in this analysis also makes the higher mode 
decomposition suspect. 



Mode#l Mode #2 Mode #3 Mode#4 Mode #5 

Op 1 1 i l— i • i I i I i I | l 1 ' 
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Mode #6 Mode #7 Mode #8 Mode #9 Mode #10 




K=3.7934 K=3.7886 K=3.7838 K-3.7779 K=3.7727 



Figure 4.6. The First 10 Normal Modes (f = 900 Hz). The Bottom Depth is 89.9 m. 
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Figure 4.7. Modal Amplitudes vs. Reduced Time - Seen. 1 1 (Range 460 m. Depth 45.7 m). 
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Figure 4.8. Modal Amplitudes vs. Reduced Time-Seen. 16 (Range 2194 m, Depth 27.4 m). 
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Figure 4.9. Modal Amplitudes vs. Reduced Time-Scen.9 (Range 4650 m. Depth 45.7 m). 
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D. SPATIAL COHERENCE (CORRELATION) 



In order to analyze the spatial coherence of the signal’s vertical structure at the VLA, 
single records from each of the three scenarios were considered. We computed the peak 
correlation of each phone’s data with a reference phone (hyd #1, closer to the bottom; 16, 
mid-array; and 32, closer to the surface) normalized such that the autocorrelation of the 
reference phone data was unity. The results provided information about how unique the 
signal was across the space (vertical direction). We can observe in Figure 4.12, from 
scenario 9 (range -5000 m), that the signal correlation decreases to values around 0.5 over 
about 5 meters in depth for the reference hydrophones 1 (bottom) and 32 (top), which 
corresponds to one element spacing and approximately two to four wavelengths. The same 
effect can be verified in Figure 4.11 which displays data from scenario 16 (range -2500 
m). These results suggest that distinct multipath arrivals at any given time are interacting 
at different locations over the array length. For scenario 1 1 (range -500 m) runs, however, 
the cross-correlation is not as low for the same difference in depth, as shown in Figure 
4.10. This is a consequence of the short distance between source and receiver which does 
not allow large variations in the signal nor many multipath arrivals. Also, by observing 
Figure 4.12, we notice that hydrophone 16’s signals are more correlated than the ones 
from hydrophones 1 and 32 for scenario 9. This could be a consequence of its relative 
position (mid-array), but the same effect was not verified in scenarios 1 1 and 16. 
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(a) Reference Hydrophone # 32 (Depth 28 m) 

♦ Phon*-1o-ption»p«#KcorT«islionlofScen»no11 -Mf 




(b) Reference Hydrophone # 16 (Depth 54.8 m) 



-lo-ptiofwpMkaDrrgtalinn Icf Scenmo 11 - MF 




(c) Reference Hydrophone # 1 (Depth 82 m) 



Figure 4. 10. Phone-to-Phone Peak Correlation - Scenario 11 (Range 460 m). 



38 




(a) Reference Hydrophone # 32 (Depth 28 m) 




(b) Reference Hydrophone #16 (Depth 54.8 m) 




(c) Reference Hydrophone # 1 (Depth 82 m) 

Figure 4.11. Phone-to-Phone Peak Correlation - Scenario 16 (Range 2194 m). 
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(a) Reference Hydrophone # 32 (Depth 28 m) 

Phorte-lo-ehone peak correlation tor Scenario 9 - MF 




(b) Reference Hydrophone #16 (Depth 54.8 m) 



.9 -MF 




(c) Reference Hydrophone # 1 (Depth 82 m) 

Figure 4.12. Phone-to-Phone Peak Correlation - Scenario 9 (Range 4960 m). 
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V. LOCALIZATION RESULTS 



This chapter presents a compilation of results obtained after a series of runs of the 
FACM algorithm. For every scenario, distinct receiver depths were tested by selecting the 
hydrophone data to be used. 

Regarding the environmental parameters used as inputs to the MMPE model, only the 
sound speed profile varied according to the respective scenario. We used only the 
bathymetry along the North radial, which is an approximate 89.9 m isobath with bottom 
composed mainly of medium to coarse grained sands. This type of sand has the following 

properties: compressional sound speed of -1800 m/s, sound speed gradient of 25s" 1 , shear 

speed of -250 m/s, density -2.03 g/cm 3 , and compressional and shear wave attenuation 
coefficients about 0.1 dB/m and 1.0 dB/m, respectively (Smith et al., 1998). During the 
experiment, the data was recorded at a sampling frequency of 4934.0 Hz. For processing, 
we extracted the received hydrophone response to the 0.1s bursts in files of 4210 points 
which corresponded to 0.8533 seconds. The reason for this procedure is that, in the 
MMPE model output, we have 512 discrete frequencies over a 600 Hz bandwidth. 

Therefore, A f = 1.171875 Hz, and we match T = ^ = 0.8533. After, we processed it 

again to reduce its length to 2048 spectral points over a bandwidth of 600 Hz to 1 200 Hz. 
The model output was zero-padded (both ends) to also contain 2048 points. 
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The spatial grid was, in most of the runs, 128 points in depth vs. 300 points in range. In 
general, it corresponded to a resolution of 0.78 m in depth and 10 m in range. For all cases, 
the source/receiver relative position proved to be a very important factor to the quality of 
the results. This can be observed in Figures 5.1 (a) and (b), where the ambiguity surfaces 
for scenario 9, based on hydrophone 16 (a), and hydrophone 2 (b), are shown. The source 
is located at depth 45.7 m, range 4650 m, and its location is depicted as a white dot. While 
the algorithm with hydrophone 16 (depth = 54.8 m) data was able to localize it within the 
some tolerance, the run with hydrophone 2 (depth = 76.3 m) data did not provide 
localization at all. 

By processing the matched filtered signal, we simulate the (idealistic) situation of a 
coherent source, and it improved the localization results a little. However, even when 
processing the raw data, we obtained satisfactory results where the main differences were 
the bigger footprint and the smaller values of the ambiguity surface. These differences can 
be observed in Figures 5.2, 5.3, and 5.4, which display results for scenario 9, scenario 16, 
and scenario 1 1 , respectively. 

Observing Figures 5.3(b) and 5.4(b), where the ambiguity surfaces from runs with MF 
data from scenarios 16 and 11 are shown, respectively, we verify that the FACM 
localization algorithm performed reasonably well (errors < 11% in range, and < 40% in 
depth). It is important to note that the sound speed profiles used were gathered four hours 
apart from experiments of scenario 9, more than two hours apart from scenario ll’s, and 
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more than 10 hours apart from scenario 16’s. The mismatch between the actual and the 
modeled environmental parameters contained in the SSP is likely significant. 

We also observed that the localization error decreased with increasing range. This was 
expected since the autocorrelation techniques are strongly dependent on the uniqueness of 
the multipath structure. The error obtained for scenario 9 (range -5000 m) tains were in the 
order of 2% in range and 30% in depth. However, scenario 11 (range -500 m) runs 
presented errors as large as 20% in range and 40% in depth. 
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(b) 

Figure 5.1. Ambiguity Surfaces for Scenario 9-(a) with Raw Data from Hyd.16 (Depth 
54.8 m), and (b) with Raw Data from Hyd.2 (Depth 76.3 m). The White Dot Depicts the 

Actual Source Location. 
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Figure 5.2. Ambiguity Surfaces for Scenario 9 - (a) with Raw Data and (b) with MF Data 
from Hyd.16 (Depth 54.8 m).The White Dot Depicts the Actual Source Location. 
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Ambiguity Surface - Raw data, Hyd.fMS, Scenario 16 




Figure 5.3. Ambiguity Surfaces for Scenario 16- (a) with Raw Data and (b) with MF Data, 
both from Hyd. 16 (Depth 54.8 m).The White Dot Depicts the Actual Source Location. 
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Figure 5.4. Ambiguity Surfaces for Scenario 11- (a) with Raw Data and (b) with MF Data, 
both from Hyd. 16 (Depth 54.8 m).The White Dot Depicts the Actual Source Location. 
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VI. CONCLUSIONS 



A. SUMMARY OF FINDINGS 

Several techniques helped us to assess the properties of the acoustic field in which we 
were applying the localization algorithm and infer their influences on it. Through spatial 
coherence analysis, we were able to verify that more unique signals led to better 
localization results, therefore, ratifying Brune (1998). Analyzing the beamformer outputs, 
we can affirm that more stratified, discrete structures, corresponding ' to the larger 
distances between the source and receiver also provided better results. Observing the 
sound speed profiles, we identified a sound channel spanning from a depth of about 30 m 
to the bottom. With the modal decomposition, we verified that the lower modes, which 
were closer to the channel axis (~ 45 m), carried more energy. When selecting a receiver 
above this channel, we observed poor results without any match, and the best results 
corresponded to hydrophones positioned closer to the sound channel axis. Another 
important fact to note is that the depths where we obtained localizations corresponded to 
the ones containing larger amplitude modes, mainly modes two to six depending on the 
source frequency chosen to be sampled. 

Regarding the use of matched-filtered or raw data as inputs to the FACM algorithm, 
our work suggests that it is not a major issue - one can have localization results applying 
raw data as useful as the ones from MF data. However, these results were based solely on 
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LFM chirp transmissions and may not apply to all signals generally. Further studies are 
needed to examine this issue more thoroughly. 



B. RECOMMENDATIONS FOR FUTURE WORK 

Once we have data from 30 of the 32 hydrophones in the VLA, an almost natural 
consideration would be the multiple-phone analysis of the localization algorithms. While 
computational load would certainly be an issue, the comparison with single-phone 
analysis, and possibly better results, would be worth it. In this case, one should consider 
the use of the array tilt information contained in the data files to correct the arrival angle, 
vertical structure, hydrophone depth, etc. 

Another possible approach would be the use of non-traditional processes, such as 
artificial neural networks (ANN) and minimum-variance spectrum estimations (MVSE), 
both in the arrival structure analysis and in the data preprocessing. This could surpass the 
over-simplifications of the traditionally used techniques and lead to better results (Ma, 
Y.L., 1997). Also, since the shallow water environment is generally a dispersive channel, 
and presents effects such as strong bottom reverberation, signal distortion and noise 
fluctuation, an analysis of the signal mismatch influences could introduce improvements 
to the process. 
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Throughout this work, only three scenarios were considered. As the NUWC-TSP 
experiment provided at least 14 other sets of data where parameters such as source 
location, signal bandwidth, and/or pulse duration differed from the ones we studied, it 
could be useful to process them to try to confirm our results. Moreover, it would allow a 
comparison of the FACM performance for longer pulses and higher frequencies. 
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